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Summary 

A pulse detonation engine (PDE) uses a series of high frequency intermittent detonation 
tubes to generate thrust. The process of filling the detonation tube with fuel and air for 
each cycle may yield non-uniform mixtures. Lack of mixture uniformity is commonly 
ignored when calculating detonation tube thrust performance. In this study, detonation 
cycles featuring idealized non-uniform H 2 /air mixtures were analyzed using the SPARK 
two-dimensional Navier- Stokes CFD code with a 7-step H 2 /air reaction mechanism. 
Mixture non-uniformities examined included axial equivalence ratio gradients, transverse 
equivalence ratio gradients, and partially fueled tubes. Three different average test 
section equivalence ratios ( <E> ), stoichiometric (0 = 1 .00) , fuel lean ( O = 0.90) , and fuel 
rich (0 = 1.10), were studied. All mixtures were detonable throughout the detonation 
tube. It was found that various mixtures representing the same test section equivalence 
ratio had specific impulses within 1% of each other, indicating that good fuel/air mixing 
is not a prerequisite for optimal detonation tube performance. 


Nomenclature 

b; Fluid element body forces, i th dimension 
c Speed of sound 

d; Dimensional characteristic boundary derivatives 

fi Mass fraction of the i th species 

h; Enthalpy of the i th species 

i re f Centerpoint of X-grid spline 

k Thermal conductivity 

kb Backward reaction rate 

kf Forward reaction rate 

n Arrhenius expression temperature exponent 

ns Number of chemical species 

p Static pressure 

p L Local tube exit pressure 

p 0 Ambient (atmospheric) pressure 

q Fluid element heat flux 

t Time 
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tcycie Time to complete one detonation tube cycle 
u X - velocity 

u; Velocity in the i th direction 

Uj Velocity in the j th direction 

ul Local tube exit velocity 

v Y - velocity 

w; Mass production of the i th species 

x Unit vector in the X direction 

x; Spatial dimension in the i th direction 

xj Spatial dimension in the j th direction 

y Unit vector in the Y direction 

Ai-D; Laminar viscosity coefficients 
A Arrhenius expression pre-exponential factor 

C p Specific heat at constant pressure 

D Binary diffusion coefficient 

D t , Thermal diffusion coefficient of the i th species 
E Total energy of flow element 

Ea Arrhenius expression activation exponent 

F Vector of convective terms from governing equations 

F Vector of convective terms from governing equations, transformed coordinates 

G Vector of diffusive terms from governing equations 

G Vector of diffusive terms from governing equations, transformed coordinates 

H Detonation tube height 

H Vector of source terms from governing equations 

H Vector of source terms from governing equations, transformed coordinates 

I Total impulse 

I sp Specific impulse 

J Jacobian Matrix 

Ke q Equilibrium species concentration 

L Detonation tube length 

L; Dimensional characteristic boundary wave amplitude terms 
M; Molecular weight of the i th species 
Mj Molecular weight of the j th species 

P CJ Chapman-Jouget point pressure 

P 0 Detonation tube pressure prior to combustion 
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Pp Plateau pressure 

Pr Prandtl Number 

P vn Von Neumann spike pressure 

R u Universal gas constant 

Sc Schmidt number 

T Static temperature 

U Vector of dependant variable from governing equations 

U Vector of dependant variable from governing equations, transformed coordinates 

V Velocity vector 

V CJ Chapman-Jouget equilibrium detonation velocity 
Vi Diffusion velocity of the i th species 

V j Diffusion velocity of the j th species 

Wfuei Detonation tube fuel weight per cycle 
X Length dimension 

Xj Mole fraction of the i th species or grid-point X-coordinate 
Xj Mole fraction of the j th species 

Y Height dimension 

8y Dirac delta function 

A, Dimensional characteristic boundary wave velocities 

/I Mixture viscosity 

p ; Species viscosity of the i th species 

/I j Species viscosity of the j th species 

(p Local equivalence ratio 

(Py Mixture viscosity weighting factor 

ri 2 nd transformed grid coordinate 

p Density 

p, Species density of the i th species 

pj Species density of the j th species 

t, 1 st transformed grid coordinate 

Ty Fluid element surface stress ( = r) 

O Detonation tube test section average equivalence ratio (Alternately given as “Phi” 
in various numerical results plots) 
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1.0 Introduction 


1 . 1 Motivation of Study 

The concept of a pulse detonation engine (PDE) for aerospace propulsion system 
applications is not new. The work at the University of Michigan in the 1950’s (Nichols 
et al., 1957) is a prime example. However, it was not until repetitive detonations with 
gaseous hydrocarbon fuels at frequencies at or above 25 detonations per second were 
demonstrated at the Naval Post-Graduate School in the mid-1980’s that it became 
apparent that practical devices might be possible, as described by Coleman (2001). 

These advances and others over the past decade have stimulated interest in the 
mainstream propulsion development community (Kailasanath, 1999; Kaemming, 2001). 
Both government and industry entities are currently considering the application of PDE’ s 
and derivative configurations to missions as diverse as Single-Stage-To-Orbit (SSTO) 
launch vehicles, rocket upper stages, commercial aircraft, missiles, and tactical aircraft 
(Coleman, 2001). At the heart of each of these unique propulsion systems is a series of 
high frequency intermittent (pulse) detonation tubes. 

PDE’s are attractive for several reasons. First, Kailasanath and Patniak (1999), 
among others, have shown that the thermal efficiency of the detonation cycle can 
approach two times the thermal efficiency of the standard Brayton cycle found in gas 
turbine engines, 49% versus 27% in an example case where both combustion processes 
began at the same initial pressure. Second, since static pressure increases substantially 
during detonative combustion, instead of slightly decreasing as in deflagrative 
combustion, it is possible to decrease the amount of compression required before the 
combustion process, thus leading to simplified, lighter, lower cost engine architectures. 
For some missions, such as air-launched missiles, it is possible to eliminate all 
mechanical pre-compression and operate the PDE as a “supercharged” ramjet, even at 
relatively low subsonic Mach numbers. Lastly, a PDE detonation chamber does not have 
to be round, nor do the detonation chambers have to be grouped in circular arrays if 
turbomachinery components have been eliminated, thus leading to a more 
aerodynamically efficient airframe integration than with a gas turbine engine. These 
PDE advantages will become more apparent in the following sections as the detonation 
tube architecture and engine gas dynamics are explored. 


1.2 PDE Architecture Description 

A variety of architectures have been proposed for pulse detonation engines, 
differing materially in tube geometry, air valve design, fuel injector placement, and other 
operational features. Figure 1 shows a generic configuration for an airbreathing pulse 
detonation tube. To begin a new cycle, the air valves on the left are opened to admit air 
from the inlet plenum. Fuel is injected into the incoming air stream as it enters the 
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initiator tube and the main detonation tube. Typically, the start of fuel injection is 
delayed slightly after the air valves are opened in order to allow a thin buffer of unfueled 
air to form between the exhaust products from the previous cycle and the fresh charge, 
preventing pre-ignition. After the inlet fuel and air valves are closed, combustion is 
initiated using the spark plug in the initiator tube. The ensuing deflagration quickly 
accelerates into a detonation, aided by the small diameter of the initiator tube and 
typically by some amount of oxygen enrichment in the initiator tube. 



Figure 1 - Basic Pulse Detonation Combustor Tube Architecture 


Once formed, the detonation wave travels down the initiator tube at a 
characteristic detonation wave speed in excess of 1500 m/s, depending on the fuel and 
oxidizer combination used and the initial conditions (pressure, temperature, 
stoichiometry). The detonation wave then propagates out of the initiator tube, detonating 
the main charge. Care must be taken in selecting appropriate geometric parameters (area 
ratio, rate of area change, cross section) that will allow the transition of the detonation 
wave into the larger tube. In some cases, the detonation may initially fail 
(shock/combustion decoupling) upon exiting the initiator tube, but will then re-initiate 
from hot spots in the main tube caused by intersecting shocks. 

Two additional features of the detonation tube shown in Figure 1 that aid in the 
initiation and propagation of the detonation wave are the convergent-divergent (C-D) 
nozzle and the actively cooled combustor walls. The nozzle, when employed, is used to 
provide a higher initial (pre-detonation) pressure in the tube by restricting the outflow 
from the tube during the fill process. This higher initial pressure aids in the initiation of a 
detonation, since detonability is a function of both pressure and temperature (Kuo, 1986). 
A higher initial pressure also increases thrust potential due to the corresponding increase 
in peak and plateau pressures caused by the pressure multiplying effect of the detonation 
wave. In the common configuration, the detonation cycle is also enabled by the presence 
of the actively cooled combustor walls. The wall temperature must be maintained below 
the autoignition temperature of the fuel/air mixture being used to prevent deflagrative 
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mode ignition prior to the initiation of the detonation. However, any mixing or buffering 
scheme that keeps the flammable mixture away from the wall could also be used to 
prevent autoignition. 


1.3 Detonation Wave and PDE Wave Cycle Descriptions 

Once fully initiated, the detonation wave assumes the characteristic Zeldovich, 
Yon Neumann, Doring (ZND) structure, shown in the fixed, or laboratory, coordinates in 
Figure 2, taken in part from Kuo (1986). The leading shock compresses the fuel/oxidizer 
mixture (A) to a pressure and temperature sufficient for rapid combustion (B). This point 
is known as the Von Neumann spike, and is approximately equal to the pressure and 
temperature rise observed across a frozen shock moving at the same Mach number as the 
detonation wave. After a short kinetic delay, combustion occurs with a characteristic 
decrease in static pressure. Point (C) corresponds to the point at which combustion is 
complete, know as the Chapman- Jouget (CJ) point. The CJ point is uniquely defined by 
the pre-detonation conditions and represents the point at which the local Mach number in 
the combustion products equals 1 in a reference frame fixed to the detonation wave. 
While it is possible to have a detonation which does not end in the CJ point, such 
detonations are generally unstable, and so the leading shock will typically self-adjust its 
speed (and therefore strength) so as to create the CJ point at the end of the combustion 
zone. The resulting shock Mach number, referenced to the speed of sound in the 
unbumed gas, is known as the CJ Mach number, and is again uniquely defined by the 
unbumed conditions. 



Figure 2 - ZND Detonation Structure in a Tube (Laboratory Coordinates) 
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In order to maintain continuity in the closed tube, a series of expansion waves (D) 
propagate toward the closed end of the detonation tube from the detonation wave front to 
reduce the pressure from the CJ pressure to what is commonly referred to as the plateau 
pressure (E), reconciling the high velocity flow following the detonation wave with the 
zero velocity requirement at the closed end wall. The tube then exhausts through the 
open end, with the blow down characteristics governed by the reflected expansion and 
compression waves propagating through the tube. Once the tube has reached a pressure 
sufficiently close to ambient pressure, the air inlet valves are opened to begin recharging 
the detonation tube leading to the initiation of a new cycle. The overall detonation cycle 
is shown graphically in Figure 3. 




(T) Detonation Initiated 





(T)New Fuel/ Air Mixture Introduced, 
Purging Out Remaining Hot Products 


(?)Tube Blows Down To 
Ambient Pressure 





(7)Dctonation Propagates 

Through Length Of Main Tube 


Figure 3 - Pulse Detonation Combustor Operational Cycle 


1.4 Scope of Numerical Study 

It is likely that in any application of detonative combustion for propulsion, the 
fuel and air will not be uniformly mixed, as is often assumed. Non-uniformity in fuel 
distribution may be created intentionally or may result from hardware limitations, and 
may be either axial or transverse, or both, depending on the source. Possible sources of 
non-uniform fuel distribution include, but are not limited to, 

• Combustion products/fresh charge buffering (purge) 

• Engine throttling 

• Emissions control 

• Wall heat transfer control 

• Detonation wave shape/strength control 

• Air and fuel valve transients 

• Non-uniform inlet air flow in space and time 

The study that follows examines the effects of various idealized fuel distribution 
non-uniformities on detonation tube thrust performance (specific impulse). Kailasanath, 
et al. (2000) studied the performance effects of unfueled purge fraction, alternately 
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referred to as partial fill, showing an increase in I sp with increasing purge fraction 
(decreasing fill fraction). Their study was conducted using a transversely uniform 
stoichiometric mixture of ethylene and air in the fueled portion of the tube, leaving the 
exit end of the tube unfueled. No spatial transition (gradient) in equivalence ratio was 
used between the two regions. Purge fractions of up to 80% were studied. This current 
study builds upon and expands the consideration of fuel distribution effects by looking at 
the effects of equivalence ratio gradients in both the axial and transverse directions as 
well as comparable partial fill cases. The numerical results are intended to be 
comparative in nature, looking for trends and relative magnitude effects between the 
different cases. The actual values of performance are not validated by test data, and so 
may be subject to different levels of offset error. However, it is expected that the 
physics-based numerical modeling techniques employed will provide correct insight into 
the relative performance effects of the different fuel distributions. 


1.5 Paper Organization 

After briefly introducing the basic concept and operation of a PDE in the 
preceding sections, this thesis presents the specific geometry, operational conditions, and 
boundary conditions being addressed in the computational study, including all the unique 
formulations added to the existing SPARK code. The SPARK code itself is then briefly 
described, along with the basic implementation of the computational grids used in the 
study. Next, the test cases used to check the accuracy of SPARK for calculating H 2 /air 
detonations are described and the corresponding results presented. 

With the code thus validated, the results of the various cases used for the study are 
presented and compared. These results are then summarized and appropriate 
conclusions drawn. Lastly, an appendix is provided justifying the selection of the 
chemistry model chosen for this study. 


2. Problem Formulation 


2. 1 Geometry 

The geometry of a detonation tube is generally simple outside of the transition 
section between the initiator tube and the main detonation chamber. While the geometry 
of the area transition is critical to operability of the device, this study will be confined to 
the constant-area main detonation chamber. This region represents the majority of the 
volumetric capacity and therefore energy release within the device, making it the 
dominant region in terms of thrust performance. To keep the geometry as simple as 
possible, the left hand side inflow plane (main air inlet valves - Figure 1) will be 
modeled as a solid wall, simulating a transitioned detonation with the inlet valves closed. 
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Typical proof-of-concept PDE detonation tubes being used today are in the range 
of 5 to 12 cm in diameter and 120 to 180 cm in overall length (including the initiator 
tube). For this study a height of 3.0 cm and a length of 30.0 cm will be used in order to 
limit the computational time required. This basic two-dimensional geometry is shown in 
Figure 4. 
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Figure 4 - Geometry for Numerical Study 


2.2 Simplifying Assumptions 

A number of simplifying assumptions are made in this study relative to the actual 
PDE cycle, as follows. 

2.2.1 DDT 

No attempt is made to model the deflagration to detonation transition (DDT) 
process, as the phenomena under study is primarily concerned with the steady state 
detonation propagation. As shown in Figure 4, the initiator tube itself is not modeled nor 
is the area change from the initiator to the main chamber modeled. A successfully 
transitioned detonation wave is numerically initiated through the imposition of a high 
temperature, high pressure ignition zone at the left-hand (inflow) side of the domain in a 
5.0 cm initiation/stabilization region. The high pressure, high temperature ignition zone 
is set uniformly to 150 atmospheres, 4000 K, and is the full height of the tube, but only 
0.6 mm long, with a composition of 25.5% H 2 O and 74.5% N 2 by weight. The total 
energy contained in the ignition region is 335 J, which is approximately 1.3% of the heat 
release from the full 30 cm detonation tube when it is uniformly, stoichiometrically 
fueled (26.2 KJ). The remaining 4.94 cm stabilization region is maintained as a uniform, 
stoichiometric unbumed mixture for each case to aid in establishing a successful 
detonation regardless of test case mixture properties. Although this detonation may be 
initially overdriven to insure propagation, it falls back to near the equilibrium Chapman- 
Jouget condition by the end of the stabilization zone. 
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2.2.2 Partial Cycle 


Only the detonation propagation and blowdown portions of the full PDE cycle, 
portions @ and © in Figure 3, are modeled. It is assumed that the purge and refresh 
cycles are identical for each case. While this might not be rigorously true in practice, 
such differences would be expected to be of significantly less importance. It is further 
assumed that the blowdown of the tube is complete when the tube exit velocity reaches 
zero anywhere along the exit plane. In a real cycle, some backflow is likely to occur due 
to the sub-ambient pressure in the tube from the last reflected expansion wave. This 
portion of the cycle is important when modeling the filling of the tube, both in terms of 
the fill time and the performance penalty, but this issue is beyond the scope of this study. 


2.3 Boundary Conditions 

The left-hand, top and bottom walls of the two-dimensional detonation tube are 
all treated as impermeable, dictating that the pressure derivative, the species 
concentration derivatives, and the normal component of velocity be set to zero at each 
wall. No-slip wall and adiabatic conditions are set at each surface by setting the 
temperature derivative and parallel component of velocity to zero. The adiabatic wall 
condition was selected to simplify the analysis by not requiring resolution of the 
boundary layer region, however the no-slip wall was maintained to better capture any 
wave reflections that might occur. 

The tube exit boundary conditions for this type of calculation have been a source 
of considerable study, with various approaches having been used by different individuals, 
such as Kailasanath and Patnaik (1999) and Wilson and Paxson (2001). Over the course 
of each cycle, there are six segments of time that need to be addressed when considering 
the outflow boundary conditions. 

1 . During the wave propagation segment, the detonable mixture between the 
detonation wave and the tube exit is undisturbed by the approaching detonation 
wave. In all of the cases considered in this study, the tube fill pressure matches 
the ambient pressure, so there is no outflow during this time segment. In most, if 
not all, real cases, there will be some low velocity flow into or out from the tube. 

2. As the detonation wave passes out of the tube, the flow at the exit plane quickly 
accelerates to a supersonic condition, as the flow between the leading shock and 
the CJ point is characteristically supersonic. 

3. After this short period of supersonic flow, the high pressure behind the reaction 
zone chokes the flow at the tube exit (Mach number equal to 1). 

4. The flow remains choked until a region of supersonic flow generated by reflecting 
gas dynamic waves reaches the tube exit, causing the exit flow from the tube to be 
supersonic for a second time. 

5. The flow then returns to the choked flow condition until the pressure has dropped 
sufficiently for the flow to become subsonic. 

6. Subsonic outflow throughout the rest of the blowdown. 
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Figure 5 shows this sequence of events at the detonation tube centerline from the baseline 
stoichiometric H 2 /air detonation calculation (Case 1) described in detail in Section 5. 
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Figure 5 - Time History of Detonation Tube Exit Flow Mach Number 

at the Tube Centerline 


This result is consistent with the results of Wilson and Paxson (2001) from their 
1 -dimensional calculations. 

The outflow boundary conditions required for each of these situations were 
formulated as follows. 

2.3.1 Subsonic Outflow 

For subsonic outflow, the methodology of Poinsot and Lele (1992) for subsonic 
reflecting outflows for the Navier- Stokes equations was adopted, from which the 
following notation is adapted. 

First, the basic characteristic velocity terms are calculated from the existing values at 
the boundary, with the downstream direction defined as the direction of flow. 

A 1 —u — c: The characteristic velocity of the upstream traveling wave (1) 
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X 2 = u : The characteristic convection velocity 


( 2 ) 


X 3 = « : The characteristic velocity at which the y velocity component is 

advected downstream (3) 

X 5 - u + c: The characteristic velocity of the downstream traveling wave (4) 


3v du 

Next, terms associated with the variation of the wave amplitude ( pc — ) for 

dt dt 


each velocity term are calculated using the 1 -dimensional Euler equations. The inviscid 
Euler equations are the best available approximation for the viscous (Navier- Stokes) case, 
and insure compatibility between the physical boundary conditions being calculated and 
the amplitudes of waves crossing the boundary. 


k = 



The amplitude term for the characteristic upstream traveling wave 


(5) 


-Z^2 *X 2 


dp dp 


: The amplitude term for the wave associated with the 


dX dX 

\ / 

characteristic convection velocity 

L - j: The amplitude term for the wave associated with the 

characteristic velocity at which the y velocity component is advected downstream 


= 4 — 

\ dX 


+ pc 


du 

dX 


The amplitude term for the characteristic downstream traveling wave 


( 6 ) 

(7) 

( 8 ) 


The X 4 and L 4 terms are omitted as they are used in reference to the 3rd spatial 
dimension, not required for this analysis. 

Next, the following collection of derivatives at the boundary are cast in terms of 
the wave amplitude associated terms for use in simplifying the governing equations. 


d(pu) 

dX 


T l 2 +Ul,+l,) 

c 2 


d \ 


(9) 


du 1 dp 

u + — 

dX p dX 


2 pc 


( 10 ) 
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( 11 ) 


U Z/1 A 

dx 3 4 


The derivatives d 2 , associated with the energy equation, and ci 5 , associated with the 

third spatial dimension, are not required for the subsonic reflecting outflow case and are 
therefore omitted here. 


Lastly, these derivatives are used in the reformulated version of the 2-dimensional 
Navier-Stokes equations below to advance the desired quantities to the next time step at 
the tube exit, assuming constant tangential stresses at the boundary. 


dp , d(pv) . . . . . 

d +C * l+ -jy = (continuity equation) 


( 12 ) 


d(pu) , , d(puv) dx u x 

— - — - + ud , + pd-, + — = — — (X-momentum equation) 

dt 1 H 3 dY dX 


(13) 


d(pv) 


dt 


+ vd j + pd ^ + 


a(pv 2 ) | dp _ 3t 22 
dY dY dY 


(Y-momentum equation) 


(14) 


where 


*n =At 


du 

dx 


2 

3 


f du dv 


dY 


V 


and 

*22 =A* 


2 3v 


2 du dv\ 

3| v ax + ar J 


(15) 


(16) 


By first solving for the density using equation (12), the u and v velocity terms can 
then be found using equations (13) and (14), respectively. For subsonic outflow, the 
static pressure is known, and is taken to be constant at 1 atmosphere for all the cases run 
during this study. Wilson and Paxson (2001) have shown that for tubes with a 
length/diameter (L/D) ratio greater than approximately 10, this is a valid assumption. 

This is because the residual effects of the blast wave on the static pressure in the vicinity 
of the tube exit have dissipated by the time the detonation tube unchokes during the last 
part of the blowdown, provided the tube pre-detonation pressure is approximately equal 
to the ambient static pressure. Lastly, the temperature at the exit boundary is determined 
by extrapolation from interior points assuming a constant local heat flux normal to the 
tube exit boundary. 
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2.3.2 Supersonic Outflow 

The values of p, T, u, and v are linearly extrapolated from the 2 immediate 
interior points, as described by Anderson (1995). The density, p , is then calculated 
using the equation of state. 

2.3.3 Choked Outflow 

The values of p, T, p and v are handled as in the supersonic outflow case, but 
with u limited to the local speed of sound. 


2.4 Detonation Tube Performance 


There are two fundamental methods for calculating the impulse generated by the 
detonation tube, as taken from Oates (1984). The first is to integrate the pressure along 
the left hand wall over the partial cycle modeled, subtracting the frictional pressure loss 
along the two side walls, 



(h \ 

| pdY 


V« 


Jleft 


( L 


\P 


Bu 

BY 


wall 


dX 


Jtop 


f L 


du 

BY' 


\P~L a „dX 


\dt 


Jbottom J 


(17) 


The second method is to integrate the momentum flux two-dimensionally at the tube exit 
over the partial cycle, correcting for non-ideal expansion, 


I H H 

7 = J jpu L 2 dY- j(p 0 -p L )dY dt 


(18) 


Both methods, shown per unit width of detonation chamber, will be implemented in the 
computer code and the relative results will be compared. The calculated impulse can 
then be used to calculate the specific impulse of the detonation tube. 


sp 


I 

W 

n fuel 


(19) 


where W fuel is the weight of fuel in the tube for one cycle. 


NAS A/TM— 2002-2 11712 


14 



2.5 Fuel-Air Combinations 


A stoichiometrically fueled 5.0 cm detonation initiation section is used to provide 
a stabilized detonation wave to the test section, as described in Section 2.2. The test 
section equivalence ratio is varied in either the X-direction or the Y-direction from case 
to case. Three test section average equivalence ratios, alternately denoted as <J> or “phi”, 
are used, one stoichiometric, one fuel-lean, and one fuel rich. Non-stoichiometric 
equivalence ratios of 0.9 and 1 . 1 were selected, representative of moderate levels of non- 
uniformity, but still well within expected stable detonation limits throughout the 
detonation tube. Different buffers or gradients of fuel/air mixture are used to achieve the 
desired average test section equivalence ratio. Transverse gradient cases are 
representative of a tube centerline fuel injection configuration, while axial gradient cases 
are representative of a non-constant tube filling process brought on by time varying tube 
inlet and exit conditions. Partial fill cases are representative of fuel lead (H 2 buffer) or 
fuel lag (air buffer) injection timing and provide a point of comparison to the previously 
referenced work of Kailasanath, et al. (2000). Baseline uniformly fueled cases are used 
to normalize the results of the subsequent cases. The test cases are grouped as follows, 
with O referencing the test section average equivalence ratio and 0 referencing the local 
equivalence ratio within the test section. 

Stoichiometric Combustion ( <P = 1.00 ) 

Case 1 - Baseline - uniformly fueled throughout test section 

Case 2 - Linear transverse gradient - 0 = 1.10 at the centerline; 0 = 0.90 at each wall 
Case 3 - Linear axial gradient - 0 = 1.10 at the closed end (after the initial 
stabilization zone); 0 = 0.90 at the open end 

Fuel-Lean Combustion ( <E> = 0.90 ) 

Case 4 - Baseline - uniformly fueled throughout test section 

Case 5 - Linear transverse gradient - 0 = 1 .00 at the centerline; 0 = 0.80 at each wall 
Case 6 - Linear axial gradient - 0 = 1 .00 at the closed end (after the initial 
stabilization zone); 0 = 0.80 at the open end 
Case 7 - 2.50 cm air buffer at the end of the test section; 0 = 1.00 up to the buffer air 
in the test section 

Fuel-Rich Combustion ( <E> = 1 . 1 0 ) 

Case 8 - Baseline - uniformly fueled throughout test section 

Case 9 - Linear transverse gradient - 0 = 1.20 at the centerline; 0 = 1.00 at each wall 
Case 10 - Linear axial gradient - 0 = 1 .20 at the closed end (after the initial 
stabilization zone); 0 = 1.00 at the open end 
Case 1 1 - 0.07 cm H 2 buffer at the end of the test section; 0 = 1.00 up to the buffer 
H 2 in the test section 
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Figures 6 and 7 show graphically the different axial and transverse fuel distribution 
schemes to be used in the “test section” of the detonation tube, excluding the baseline 
uniformly fueled cases (Cases 1, 4, and 8) and the partial-fill cases (Cases 7 and 11). 

Test cases performed using the opposite orientation of fuel distribution non- 
uniformity (transverse distributions that were fuel lean at the centerline and fuel rich at 
the walls, and axial distributions that increased in equivalence ratio along the length of 
the tube instead of decreasing) yielded nearly identical detonation tube performance 
results as the cases shown above, and so were not included in this study. 
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3. Numerical Tool Description 

The numerical tool chosen for this study was the SPARK 2D Navier-Stokes Code 
developed at the NASA Langley Research Center. The SPARK code was developed 
primarily for the study of high speed reacting flows, particularly those found in scramjet 
engines, under the National AeroSpace Plane (NASP) program. SPARK is a versatile 
time accurate code that has been used in a number of supersonic combustion and 
detonation studies (Drummond and Mukunda, 1988; Drummond, 1988; Carpenter, 1989). 
The code has been previously described in detail by Drummond (1990), so this 
description will be limited to a general overview of the principal features salient to this 
study. 


3 . 1 Governing Equations 

SPARK uses the strong conservation form of the Navier-Stokes, energy, and species 
continuity equations as derived by Williams (1985). 
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Continuity: 


^+v-H=o 

Momentum: 

^0+v(pvv)=w ■T + pf l f,b, 

dt i= i 
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Species Continuity: 
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(23) 

(24) 
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The expression for diffusion velocity, V j , assumes that N 2 is the major species present, 

with all other species present in “trace” amounts. It can then be assumed that only the 
diffusion of each species with N 2 be considered, and further that all binary diffusion 
coefficients are equal, as described by Drummond (1990). 

Anderson (1995) has pointed out that when the governing equations are cast in the 
strong conservation form as above, there is no issue with capturing shocks in the flow 
field. Since the equations are set up in terms of fluxes instead of the primitive variables 
(j p,u,v,p,etc .), there are no discontinuities across shock waves within the governing 
equations that require separate consideration for closure. 


3.2 Detailed Chemistry 

Following the work of Cambier and Adelman (1988) and Lynch and Edelman 
(1994), a 7-step, 7-species hydrogen-air detailed chemical mechanism was chosen for this 
study instead of the more standard 18-step, 9-species mechanism of Jachimowski (1988). 
For computational efficiency, it is desirable to utilize the simplest possible mechanism 
that can with reasonable accuracy provide the correct time evolution and chemical 
composition of the detonation wave and the resultant combustion products. The 7-step 
mechanism is a simple reduction of the 18-step mechanism through the elimination of 
HO 2 and H 2 O 2 from the chemical system, leaving H 2 , O 2 , H 2 O, OH, H, O, and N 2 (inert). 
This simplification can be justified when it is considered that in the basic ZND detonation 
model, ignition occurs after the leading shock, at which point the flow temperature is 
elevated to a point where the chemistry of HO 2 and H 2 O 2 is no longer dominant. 

The Appendix provides further analysis to justify the use of the simplified chemical 
mechanism. 

Each chemical reaction is modeled using the standard Arrhenius expression for 
the forward reaction rate, 


-Eg 

k f = A*T n * e R “ T 


(28) 


For the 7-step mechanism used in this study the values of A, n, and Ea are as follows. 
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Table 1 - 7-Step H 2 -Air Mechanism 


Reaction 

Ref.# 

Reaction 

A 

n 

Ea 

(cal/g-mole) 

1 

H 2 +0 2 <=> 2 OH 

0.1700el4 

0.00 

48150 

2 

h+o 2 ^oh+o 

0.1420el5 

0.00 

16400 

3 

h 2 +oi /<= $h 2 o+h 

0.3160e08 

1.80 

3030 

4 

0 + H 2 <=> OH + H 

0.2070el5 

0.00 

13750 

5 

OH + OH H 2 0 + 0 

0.5500el4 

0.00 

7000 

6 

OH + H + M ^ H 2 0 + M 

0.2210e23 

-2.00 

0 

7 

H + H + M <=> H 2 +M 

0.6530el8 

-1.00 

0 


The 3 rd body efficiencies of all species for reactions 6 and 7 are assumed to be one. The 
backward rate constant is then found from the simple expression, 


K 



(29) 


where K eq is found using a Gibbs free energy minimization scheme similar to those found 
in chemical equilibrium codes such as the NASA Glenn Research Center’s CEA code 
developed by Gordon and McBride (1994). The chemical system involving the 6 
reactive species is solved at each time step, thus providing the source term for the species 
equation (27). The mass of nitrogen present is determined by summing the mass 
fractions of each of the other species and then subtracting the result from the necessary 
total of 1 . Mass conservation is maintained throughout the computational domain at each 
time step. 


3.3 Thermodynamic and Transport Properties 

SPARK utilizes previously developed analytical expressions for the required 
thermodynamic and transport properties used throughout the code, which, as with the 
previously discussed aspects of the code, have been related in detail in other publications, 
such as Drummond (1990). Species specific heat at constant pressure and species Gibbs 
energy are calculated using 4 th and 5 th order polynomial curve fits of temperature, 
respectively. The species laminar viscosity, , is calculated using an expression known 
as Sutherland’s law, 
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3.4 Computational Grid 

Looking back at the detonation structure shown in Figure 2, it is easily seen that 
the region in the immediate vicinity of the detonation wave front will determine the 
required physical grid spacing for the computational model. If the goal of this study were 
to examine the detailed structure of the shock wave coupled to the reaction zone, it would 
be necessary to provide a computational grid spacing on the order of the width of the 
leading shock wave (shown by Shapiro (1953) to be on the order of 1.4x10 7 m), and 
probably finer. Fortunately, it is not necessary to resolve this structure in order to obtain 
accurate flow properties required to calculate the performance values that are the goal of 
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this study. It is sufficient to provide enough grid resolution to capture the shortest 
reaction times of the chemical system as the flow behind the shock is convected 
downstream. For the 7-step chemical system and conditions investigated in this study, 
this worked out to be on the order of 1x10 5 m. Grid spacing significantly larger than 
this, on the order of 1x10^ m , resulted in oscillations in the combustion zone. When 
even courser grids were tried, the calculation went completely unstable. From a practical 
standpoint, it is desirable to use the largest acceptable grid spacing in order to reduce the 
number of required grid points. Also, since the system of equations being solved are 
numerically “stiff’, Radhakrishnan (1991) has shown that going to a finer grid spacing 
results in a proportionally smaller required time step, leading to very long computational 
times. 


Since an X-grid spacing of 1x10 5 m would yield 30,000 grid points in the x- 
direction for each Y-grid location for a 30 cm tube, it becomes immediately apparent that 
this grid spacing needs to be confined to the region where it is required, namely the 
leading section of the detonation wave where steep gradients in all the principal flow and 
composition variables are present. Since the detonation wave is moving in space, the 
fine grid must move with it. This was accomplished for this study by starting with a 
uniform 300-point grid in X at the fine grid spacing (lxlO -5 m). As the detonation wave 
approaches the right-hand edge of the computational domain, the grid is stretched by 
approximately 30 times the fine grid spacing at the left hand side where the reaction has 
already occurred, lengthening the overall domain in X. The fine grid section and the 
coarser grid section are connected through a spline fit, 


y,.=y m +ay ci 


(ay -ay, ) 

V coarse fine / 

1 + ^-0.05555 (i-i ref ) 5 


(35) 


where i ref is the center grid point of the spline, which is set to 60 for all the simulations in 

this study. This value of i ref gives approximately 100 grid points at the fine grid spacing 

at the right hand side of the domain. The wave structure is thus always maintained within 
the fine grid section, which is effectively shifted to the right every time the detonation 
wave travels 30 times the fine grid spacing. The current result is then interpolated onto 
the new grid and the calculation is continued. Once the desired length of the full 
detonation tube is reached, the fine grid is maintained at the right hand (exit) end of the 
tube. Figure 8 shows the physical X-grid spacing at detonation wave positions (domain 
lengths) of 3 cm, 10 cm, and 30 cm. 
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X Grid Point 

Figure 8 - Variable X-Grid Spacing 


The Y-grid spacing remains the same throughout each calculation. Grid points 
are clustered near the wall in order to capture with reasonable fidelity the severe gradient 
of velocity generated by the detonation wave in that region. As with the X-grid, 
significantly finer grid spacing yields proportionally smaller time steps, leading to 
excessively long computational times. Grid spacings as small as 5x10 6 m did not 
require time steps smaller than those required for the finite-rate chemistry. The general 
Y-grid distribution used throughout the study is shown in Figure 9. Minor adjustments to 
this grid were used case-to-case to minimize run times while maintaining detonation 
wave stability. 



Y Grid Point 
Figure 9 - Y-Grid Spacing 
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Once the grid has been generated, the method of Smith and Weigel (1980), 
referred to as the “two-boundary technique,” is used to transform the mesh into a uniform 
2-dimensional boundary- fitted rectangular grid (<^,r/) for computation. The disconnected 
domain boundaries are defined by analytic functions, which are then connected together 
using parametric polynomials to generate the grid. This transformation allows for more 
efficient computation of the time-accurate solution, which will be further elaborated on in 
the next section. The solution is then transformed back into the original physical space 
for output. 


3.5 Calculation Implementation and Numerical Solution 

Once the geometric domain has been discretized into the computational grid, the 
governing equations are vectorized for solution in the form, 


dU dF dG — 

+ + = H 

dt dX dY 


(36) 


where the U vector represents the dependant variables, the F and G vectors the 

convective and diffusive terms in X and Y, respectively, and the H vector the source 
terms. In the transformed coordinate system, equation (36) becomes, 


dU dF dG Tv 
dt dt, dri 


where 


U = JU , 

07 ] 07 ] 

V_ar- ar =, 


H = JH 
and where 


(37) 

(38) 

(39) 

(40) 

(41) 
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(42) 


“ d£ dri dri 

This vectorized equation is then solved using a MacCormack explicit predictor- 
corrector scheme as described by Drummond et al. (1991), that is 4 th order accurate in 
space and 2 nd order accurate in time. While the code allows for both partially and fully 
implicit schemes to be used, these options were not selected for this study in order to 
maintain the greatest accuracy possible. The time step was limited so as to allow a 
maximum change in any species mass fraction of 0. 8% per iteration at any grid point. 
With these options selected the typical simulation contained in this study required 60 
hours to complete on a Silicon Graphics Octane2 workstation with a 400 MHz R12000A 
processor. 

4. Code Validation 

4. 1 Introduction 

Before beginning to utilize the SPARK code for the configuration and conditions 
of interest, it is necessary to demonstrate that the code can accurately compute one or 
more known solutions in similar flow regimes. The SPARK code has been previously 
validated for such cases as 

• Non-reacting 

o Laminar spatially developing shear layer 
o Turbulent spatially developing mixing layer 
o Mach 14 laminar 15 degree compression comer 
o Transverse jet behind a rearward facing step 

• Reacting 

o H 2 -air turbulent mixing and reaction 

o Burro ws-Kurkov Mach 2.44 H 2 fueled vitiated scramjet combustor 
All of these cases are described in Drummond (1991). 

Additionally, prior to initiating the detonation computations included in this 
study, several standard non-reacting test cases were performed, including a supersonic 
expansion comer, a supersonic compression comer, and a H 2 -N 2 shock tube problem. In 
all cases the results were consistent with the corresponding classical analytical solution. 
In addition to the above validations, it was necessary to perform validation calculations 
specifically for detonation. Two cases were selected for which known solutions with 
detailed kinetic mechanisms exist or could be easily calculated using other validated 
codes, as follows. 
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4.2 H 2 -Air Oblique Detonation Wave 


Thaker and Chelliah (1997) used the SPARK code to calculate the detonation 
structure of an oblique detonation wave created from a Mach 9.3 stoichiometric H 2 -air 
mixture flowing over a 27 degree wedge. The geometry and computational domain are 
shown in Figure 10. 



Figure 10 - Oblique Detonation Geometry and Computational Domain 


As shown, the computational domain was tilted 27 degrees to coincide with the 
wedge surface. A 200 by 150 uniform grid was used to model the 5 mm by 2.5 mm 
region. 0.5 mm of open boundary was included in the lower boundary before the leading 
edge of the wedge. A slip-wall, adiabatic boundary condition was used along the wedge 
wall. 


The flow entering the domain was at the following conditions, 

Pressure = 1 atmosphere Temperature = 300 K 

X-direction velocity = 3430 m/s Y-direction velocity = -1748 m/s 

H 2 mass fraction = 0.0283 0 2 mass fraction = 0.2265 

N 2 mass fraction = 0.7452 

The baseline calculation of Thaker and Chelliah (1997) was first repeated using 
the 12-step H 2 -air mechanism used in their study. This mechanism was taken from the 
full 18-step mechanism given in the Appendix by removing the H 2 0 2 reactions. The 
resulting temperature and H 2 0 mass fraction results are shown in Figure 1 1 and 
Figure 12. These results are consistent with those of Thaker and Chelliah (1997). 
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Figure 12 - 12-Step Mechanism Oblique Detonation H 2 0 Mass Fraction Contours 


Next, the calculation was repeated using the 7-step mechanism chosen for this 
study. As shown in Figures 13 and 14, the initial shock/reaction zone is very similar to 
that obtained with the 12-step mechanism. However, it is apparent that the final stages of 
heat release are somewhat slower with the 7-step mechanism. This is consistent with the 
1-D results to be shown in the next section. It is not expected that this slight delay in heat 
release will influence the comparative results of this study. 
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4.3 Uniform H 2 -Air Mixture ZND Detonation 


As a final code validation case, the detonation structure computed by SPARK in a 
uniform stoichiometric H 2 -air mixture initially at 1 atmosphere pressure and 300 K (Case 
1 of this study, detonation tube centerline results) was compared to that obtained from the 
1 -dimensional ZND code developed by Shepherd (1986). To run the ZND code, the 
desired kinetic mechanism is first input to the CHEMKIN code, described by Kee et al. 
(1989), that provides the kinetic and thermodynamic properties required. Next, the 
Gordon and McBride (1994) CEA code, or any other equilibrium detonation code, is used 
to determine the detonation (CJ) velocity. With these inputs, ZND calculates the 1-D 
detonation structure by first assuming a normal shock based on the detonation velocity 
and then reacting the mixture via the kinetic mechanism and gas properties from 
CHEMKIN using the post-shock conditions as the initial conditions. ZND performs 
these calculations in the reference frame fixed to the detonation wave, shown in Figure 
15, so the velocity must be transformed back to the stationary reference frame by adding 
the detonation wave velocity. All of the static properties are unaffected by the change in 
reference frames. 


Unburned Gas 
Vi=0 



© 


© 


Detonation 

Front 


Burned Gas 

◄ v 2 




Stationary 

Detonation 

Front 


Unburned Gas 
Vj = V c j ► 



© 


© 


Burned Gas 
► V 2 


Stationary Reference Frame 
(SPARK) 


Reference Frame Fixed to Detonation Wave 
(ZND) 


Figure 15 - Stationary and Detonation Wave Reference Frames 

Figures 16 through 22 show the values of pressure, temperature, velocity, and 
H 2 0, OH, O, and H mass fractions calculated by SPARK against those calculated by 
ZND using both the 7-step kinetic mechanism and the full 18-step mechanism. The 18- 
step ZND results are shown to additionally quantify the accuracy of the 7-step 
approximation beyond what is shown in the Appendix. The results from ZND begin after 
the leading shock. The results shown for SPARK are taken from the baseline case, Case 
1, which will be discussed further in Section 5. As shown, the ZND and SPARK results 
agree very well. The slight delay in completing the heat release that was discussed in the 
previous section is again evident in the temperature and H 2 0 mass fraction plots. This is 
a very minor concern. Based on these results and those preceding, it appears that SPARK 
with the 7-step H 2 -air mechanism is adequate for the desired detonation calculations. 
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Distance (cm) 

Figure 17 - Comparison of Temperature Profiles Between SPARK and ZND Codes 
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Figure 19 - Comparison of H 2 O Mass Fraction Profiles Between SPARK 

and ZND Codes 
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Figure 20 - Comparison of Detonation OH Mass Fraction Profiles Between SPARK 

and ZND Codes 



Figure 21 - Comparison of Detonation O Mass Fraction Profiles Between SPARK 
and ZND Codes 
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Figure 22 - Comparison of Detonation H Mass Fraction Profiles Between SPARK 

and ZND Codes 


5. Computational Results 

Having shown the general applicability of the SPARK code with the 7-step H 2 /air 
mechanism for detonation calculations, we will now look at the results of each of the 1 1 
cases outlined in Section 2.5. Figure 4 is repeated below as Figure 23 as a reminder of 
the geometry and detonation initiation/stabilization approach used throughout the study. 


a 

o 

m 


1 
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Test Section - Tailored Mixtures 


5 cm 
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Figure 23 - Geometry for Numerical Study 
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5.1 Introduction 


Before moving on to compare the results of the subject test cases, there are a 
number of introductory questions that need to be addressed. These questions deal with 
the stability of both the numerical solution and the detonation itself, the appropriate 
impulse calculation to use in the I sp calculation, and the effect of equivalence ratio on I sp 
for uniform mixtures. It will also be useful to look at simplifications that can be made in 
the analysis of the 2-dimensional numerical results. These issues can all be addressed by 
looking at elements of the solutions that will be more fully presented in later sections. 

5.1.1 Dimensionality of Results 

It is appropriate to focus on the last issue first, namely, any simplifications that 
can be made in analyzing the 2-dimensional results obtained from SPARK. Detonation 
tubes are often simulated using 1 -dimensional codes, as demonstrated by Kailasanath and 
Patniak (1999), Wilson and Paxson (2001), and Sterling et al. (1995). If the assumption 
of primarily 1 -dimensional behavior for uniform and axially varying mixtures can be 
verified, then those cases can be visualized and compared by looking at 1 -dimensional 
slices along the length of the tube, typically taken along the centerline. It is obvious that 
the transversely varying mixtures will still need to be examined to some extent using 2- 
dimensional methods, but even these may be looked at 1 -dimensionally for some flow 
variables. 

Figures 24 to 28 show 2-dimensional plots of pressure, transverse velocity and 
temperature from Cases 1, 2, 3, 4, and 7 as examples of a completely uniform mixture, a 
transversely varying mixture, an axially varying mixture, a uniform mixture with a step 
change in equivalence ratio along the tube, and a partially fdled tube, respectively. All 
the plots in Figures 24 to 28 represent the flowfield just prior to the detonation wave 
exiting the tube. While all the cases show small traces of transverse velocity, only Case 
2 shows any signs of significant multi-dimensionality, seen in the temperature plot for 
that case. However, the pressure plot for Case 2 is still very 1 -dimensional, which is to 
be expected, since there can be no net flow in the transverse direction. Thus for all the 
cases, one-dimensional comparisons of pressure should be valid. 

It is of interest that there is no observed boundary layer in any of the plots. A 
boundary layer is present in all cases, but it is so thin as to be unobservable in the scale of 
the plots, and it has essentially no effect on the flow features of the cases under study. It 
should also be noted that in Case 7, a partial fill case, the detonation wave has propagated 
into the air only region at the end of the tube, leading to the region of lower temperature 
at the end of the tube seen in that plot. 
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Figure 24 - Case 1 Pressure, Transverse Velocity, and Temperature Contours 
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Figure 25 - Case 2 Pressure, Transverse Velocity, and Temperature Contours 
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Figure 26 - Case 3 Pressure, Transverse Velocity, and Temperature Contours 
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Figure 27 - Case 4 Pressure, Transverse Velocity, and Temperature Contours 
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Figure 28 - Case 7 Pressure, Transverse Velocity, and Temperature Contours 


5.1.2 Numerical and Detonation Stability 

It is of great importance to establish that each simulation is both numerically and 
physically stable. This is most easily accomplished by looking at the detonation wave 
speed. The wave speed was determined by tracking the wave position, defined as the 
tube centerline position with the steepest X-velocity gradient, as a function of time during 
the wave propagation period of the calculation. The wave speed was found to be very 
sensitive to both concerns during the validation stages of this study. The equilibrium 
value of the Chapman- Jouget wave speed is easily calculable using the CEA code of 
Gordon and Mcbride (1994). In Figure 29, the three baseline case wave speeds are 
compared to each other, from which it can be seen that each simulation is stable, and that 
wave speed is proportional to equivalence ratio. Figures 30 through 33 show the wave 
speed for each baseline case compared to several reference detonation wave speeds 
calculated using the CEA code. After being initially overdriven to insure detonation 
initiation, the wave relaxes to just under the CJ wave speed in each case, and then slowly 
converges on the equilibrium CJ wave speed throughout the “steady state” propagation. 
The steady state wave speed is generally within 2% to 3% of the equilibrium value, also 
indicative of a stable detonation. 
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Figure 29 - Detonation Wave Speeds for Cases 1, 4, and 8 
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Figure 30 - Case 1 (0 = 1.0) Wave Speed Compared to Equilibrium CJ Wave Speeds 
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Time (s) 

Figure 31 - Case 4 (O = 0.9 ) Wave Speed Compared to Equilibrium CJ Wave Speeds 



Time (s) 


Figure 32 - Case 8 ( <E> = 1 . 1 ) Wave Speed Compared to Equilibrium CJ Wave Speeds 
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5.1.3 Total Impulse Calculation 

As discussed in Section 2.4, there are two methods for calculating the impulse 
generated by the detonation and subsequent blowdown. The first is to integrate the 
pressure acting on the head-end of the tube over the cycle time, subtracting off the 
viscous drag of the walls, as given in equation 17. It is also acceptable to integrate the 
momentum flux exiting the tube over the cycle time, correcting for the incomplete 
expansion to ambient pressure, as given in equation 18. Figure 33 shows the results of 
both impulse calculations versus time for Case 1 . It is not possible to differentiate the 
impulse generated by the initiation and stabilization zones, so calculated total impulse 
includes the effects of these regions. The weight of H 2 in the stabilization zone is 
included in the total fuel weight in order to calculate the overall specific impulse. 
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Figure 33 - Impulse Calculations for Case 1 - Uniform Stoichiometric Mixture 

From Figure 33 it is seen that the head end pressure calculation has come close to 
its maximum value, while the exit momentum flux calculation is still rising. The two 
impulse values are about 6% apart at the termination of the calculation due to flow 
reversal at the tube exit boundary. This result is understandable in that flow is still 
exiting the tube at a considerable velocity at this point (see Figure 5). However, the 
pressure at the head-end has returned close to the initial pressure. The head-end pressure 
will in fact go below the ambient pressure if the calculation is allowed to continue due to 
reflected expansion waves in the tube. It is expected that if the calculation were run out 
far enough so that the tube returned to a quiescent state, these two impulse values would 
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equalize. However, the point at which the head-end pressure begins to go sub- 
atmospheric is precisely the point at which it would be expected that the inlet air valves 
would open, drawing in fresh fuel and air and preventing the head pressure from going 
significantly negative. Since the thrust wall pressure integration has approached a 
maximum, it is evident that this is the value to use in the I sp calculation, being the closest 
to the expected maximum value. All further references to I sp values are based on the 
head-end pressure integration method of equation 17. 

Before leaving Figure 33, it is interesting to note that the wall shear force impulse 
penalty also included with the figure shows that the wall drag only represents a penalty of 
approximately 0.75% of the total impulse. While the viscous drag is included in the total 
impulse calculation, it could be readily ignored without significantly altering the results 
of any of the study cases. 

5.1.4 Comparison of Baseline Uniformly Fueled Cases 

It will be useful in understanding the effect of equivalence ratio gradients if the 
general effect of equivalence ratio in uniformly fueled mixtures in the range of interest is 
known. Figure 34 shows that I sp is inversely proportional to equivalence ratio, with the 
lean case achieving the highest I sp (4706 s), the stoichiometric case in the middle (4450 
s), and the fuel rich case yielding the lowest I sp (4177 s). If these maximum I sp values 
are plotted against overall equivalence ratio, as in Figure 35, it is seen that the I sp 
dependence on equivalence ratio in this range is approximately linear. Inclusion of the 
stabilization zone raises the overall equivalence ratio of the fuel lean cases to 0.92, and 
lowers the overall equivalence ratio of the fuel rich cases to 1 .08. Since the effect of the 
stabilization zone cannot be isolated within the calculation, this equivalence ratio bias is 
included in all Isp calculations. 

Since I sp is calculated from the force generated by the fluid pressure pushing 
against the closed end of the tube, the head-end pressure can be plotted against time to 
get more insight into the source of the variation in I sp with equivalence ratio, as in 
Figure 36. It is evident that decreasing the equivalence ratio below an equivalence ratio 
of one decreases the plateau pressure achieved, but not in proportion to the decrease in 
the amount of fuel used, thus causing an increase in I sp for Case 4. However, no 
significant increase in plateau pressure was achieved by over- fueling the tube as in Case 
8, causing a decrease in I sp due to the increase in fuel used. It is of interest to note that 
the duration of the pressure plateau and the slope of the blowdown did not change 
significantly as a function of equivalence ratio. 
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Figure 34 - I sp for Baseline Uniformly Fueled Cases 1, 4, and 8 



Figure 35 - I sp as a Function of Overall Equivalence Ratio, including Stabilization Zone, 

for Cases 1, 4, and 8 
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Figure 36 - Detonation Tube Head-End Pressure for Cases 1, 4, and 8 

Having dealt with these introductory issues, we can now move on to the analysis 
of the performance results for each group of mixtures. 


5.2 Stoichiometric Mixtures 

Before looking at the performance results, it is of value to verify the numerical 
and physical stability of the detonations. Figures 37 through 40 show the wave speed 
results for the first 3 cases, first compared to each other and then to several reference 
equilibrium values calculated using the CEA code. Each simulation is clearly stable. It 
is of interest to note that the transverse fuel distribution, Case 2, reaches the same wave 
speed as the uniformly fueled mixture, Case 1, after initial regions of overshoot and 
undershoot, and that this wave speed is stable. This indicates that the detonation wave is 
underdriven in the rich mixture and overdriven in the lean mixture. We will come back 
to this point after looking at the performance results. The wave speed of the axially 
varying fuel distribution case follows the initial equivalence ratio distribution fairly well, 
deviating from the equilibrium value slightly more at the step change in equivalence ratio 
located at the 5 cm point (~2.2xl 0 -5 s) as there is not enough time for the detonation 
wave to react to the change in mixture at that point. 
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Figure 37 - Detonation Wave Speeds for Stoichiometric Cases 1, 2, and 3 



Figure 38 - Case 1 (Uniform Mixture) Wave Speed Compared to 
Equilibrium CJ Wave Speeds 
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Figure 39 - Case 2 (Transverse Gradient) Wave Speed Compared to 
Equilibrium CJ Wave Speeds 



Time (s) 

Figure 40 - Case 3 (Axial Gradient) Wave Speed Compared to 
Equilibrium CJ Wave Speeds 
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Figure 41 shows the specific impulse traces versus time for each of the 
stoichiometrically fueled cases. In the scale of Figure 41, there is essentially no 
difference in I sp across the three simulations. Figure 42 provides a closer view of the 
final section of the I sp curves from Figure 41. We can see that there is some difference in 
the final values of I sp , but that the results are all within 1% of each other, with the axial 
distribution slightly outperforming the other 2 distributions. 

Figure 43 shows that the head-end pressure is essentially unchanged for the three 
stoichiometric cases. Case 3 shows a slight increase in the head end pressure during the 
plateau region, consistent with the slightly improved performance shown in Figure 42. 
The plateau duration and blowdown contour are virtually identical for all three cases. 



Figure 41 - Specific Impulse of Stoichiometric Cases 1, 2, and 3 
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Figure 42 - Enlarged View of Final Specific Impulse of Stoichiometric 

Cases 1, 2, and 3 



Figure 43 - Head-End Pressure as a Function of Time for Stoichiometric 

Cases 1, 2, and 3 
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Lastly, we can look at plots of H 2 or H 2 0 mass fraction to determine if any 
significant mixing or post-detonation combustion is occurring in the detonation tube. 
Figure 44 shows the hydrogen mass fraction in the tube at five different times during the 
simulation. The first three plots are during the detonation propagation, and the last two 
plots are during the blowdown. The post-detonation H 2 distribution remains essentially 
stratified throughout the entire propagation and blowdown, with the composition 
distribution at the left of the tube stretching to the right to fill the tube as the blowdown 
progresses. The products of the uniformly fueled stoichiometric stabilization zone 
eventually push out most of the stratified combustion products by the end of the 
simulation. The relatively high levels of residual H 2 at the left-hand (closed) end of the 
tube seen in this time sequence are a result of the thermal breakdown of H 2 0 from in the 
high pressure, high temperature region used to initiate the detonation and are not a 
product of the detonation process itself. It can be assumed that negligible mixing occurs 
in the axial direction if little mixing is observed in the transverse direction, since the axial 
composition gradients are approximately 1/1 6 th as steep as the transverse composition 
gradients. This is bom out by the one-dimensional time history of H 2 0 mass fraction for 
Case 3 shown in Figure 45. 



Figure 44 - Two-Dimensional H 2 Mass Fraction Distributions for Case 2 - 

Transverse Gradient 


NAS A/TM— 2002-2 11712 


48 







Figure 45 - One-Dimensional H 2 O Mass Fraction Distributions for Case 3 - 

Axial Gradient 


5.3 Fuel Lean Mixtures 

The wave speed results for the fuel lean cases, Figures 46 through 50, are much the 
same as the stoichiometric cases. The uniformly fueled mixture, Case 4, and the 
transversely varying mixture, Case 5, follow a similar profde of wave speed. The axially 
varying mixture, Case 6, follows the initial mixture profile at near the equilibrium 
detonation speed. The added case, a uniform stoichiometric mixture with an air buffer at 
the end of the tube, shows the expected drop off in wave speed as the detonation fails in 
the unfueled air. The wave speed falls off precipitously as the detonation wave enters the 
unfueled region, limiting the transfer of energy to the buffer air region. 
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Figure 46 - Detonation Wave Speeds of Fuel Lean Cases 4, 5, 6, and 7 
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Figure 47 - Case 4 (Uniform Mixture) Wave Speed Compared to Equilibrium 

CJ Wave Speeds 


NASA/TM — 2002-2 11712 


50 



Time (s) 


Figure 48 - Case 5 (Transverse Gradient) Wave Speed Compared to 
Equilibrium CJ Wave Speeds 



Time (s) 

Figure 49 - Case 6 (Axial Gradient) Wave Speed Compared to 
Equilibrium CJ Wave Speeds 
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Figure 50 - Case 7 (Partial Fill (Air)) Wave Speed Compared to Equilibrium 

CJ Wave Speeds 


The specific impulse results, Figure 51, are also similar to the stoichiometric results, 
allowing for the difference in magnitude due to the change in overall equivalence ratio as 
documented in Section 5.1. The total variation in specific impulse is on the order of 1%. 
We again see the axial distribution slightly outperforming the uniform and transverse 
distributions. The new “partial fill” case slightly under-performs the others. This 
indicates that it would be slightly better to throttle a PDE by decreasing the equivalence 
ratio through reduced fuel flow rather than by partial filling, if that is a practical approach 
when other factors are considered. 

The head end pressure time histories in Figure 52 show more case-to-case 
variation than did the stoichiometric mixtures. The axially varying mixture has a 
continually decreasing plateau pressure, and the partial fill case has a shortened pressure 
plateau. While there are differences, it is evident from the specific impulse results that 
the positive and negative aspects of these differences work to cancel each other out. 
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Figure 51 - Enlarged View of Final Specific Impulse of Fuel Lean 

Cases 4, 5, 6, and 7 
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Figure 52 - Head-End Pressure as a Function of Time for Fuel Lean Cases 4, 5, 6, and 7 
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Again checking for mixing in the detonation products for the transversely varying 
mixture, Case 5, Figure 53 shows a time series of two-dimensional plots of O 2 mass 
fraction. Since there are no fuel rich zones in this distribution, it is expected that there 
will be no more than trace amounts of H 2 present in the combustion products, so O 2 
serves as a better indicator of mixing and/or combustion within the tube. The same 
characteristic behavior is seen here as was seen for the stoichiometric case. The stratified 
post-detonation composition distribution remains stratified as it is pushed out of the 
detonation tube during the blowdown process. No significant mixing or diffusion leading 
to post-detonation combustion is observed. 
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Figure 53 - Two-Dimensional O 2 Mass Fraction Distributions for Case 5 - 

Transverse Gradient 


5.4 Fuel Rich Mixtures 


The wave speed results, Figures 54 through 58, for the fuel rich cases also show 
similar behavior to those for the previous two overall equivalence ratios. The uniformly 
fueled mixture wave speed, Case 8, and the transversely varying mixture wave speed, 
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Figure 54 - Detonation Wave Speeds of Fuel Rich Cases 8, 9, 10, and 1 1 
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Figure 55 - Case 8 (Uniform Mixture) Wave Speed Compared to Equilibrium 

CJ Wave Speeds 
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Figure 56 - Case 9 (Transverse Gradient) Wave Speed Compared to 
Equilibrium CJ Wave Speeds 



Figure 57 - Case 10 (Axial Gradient) Wave Speed Compared to 
Equilibrium CJ Wave Speeds 
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Figure 58 - Case 1 1 (Partial Fill (H 2 )) Wave Speed Compared to Equilibrium 

CJ Wave Speeds 


Case 9, do not agree quite as well as before, but are still relatively close. The axially 
varying mixture wave speed, Case 10, still follows the initial mixture distribution. The 
partial fill distribution, Case 1 1 , this time with H 2 at the end of the tube instead of air, 
shows a marked increase in wave speed as the detonation wave enters the end buffer 
region. This is due to the change in the local speed of sound that is significantly higher in 
hydrogen than in air. 

The specific impulse results (Figure 59) are, if anything, more tightly clustered 
than in either of the previous groups of mixtures. The uniformly fueled case and the 
axially varying case are almost identical, with the transversely varying case slightly 
lower. The partial fill case has gone from being the worst performer with an air buffer 
for the fuel lean cases, to being the best performer for the fuel rich cases. However, again 
the total variation in results is on the order of only 1%. 

Looking at the head end pressure time histories in Figure 60, we see almost no 
variation at all, despite the slightly increased variability in wave speed observed in 
Figure 54. 

Lastly, we see in Figure 61 the same lack of mixing in the post-detonation 
composition distribution of the 2-dimensional plots for the transverse mixture in Case 9. 
H 2 mass fraction has again been used for these plots since there are fuel rich zones in the 
pre-combustion mixture. 
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Figure 59 - Enlarged View of Final Specific Impulse of Fuel Rich 
Cases 8, 9, 10, and 11 
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Figure 60 - Head-End Pressure as a Function of Time for Fuel 
Rich Cases 8, 9, 10, and 1 1 
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Figure 61 - H 2 Mass Fraction Distributions for Case 9 - Transverse Gradient 


6. Summary and Conclusions 


6.1 Summary of Results 


A total of eleven 2-dimensional detonation tube simulations were conducted. Of 
these, three were run with stoichiometric Fl 2 /air mixtures, four were run with fuel lean 
H 2 /air mixtures, and four were run with fuel rich H 2 /air mixtures. An approximately 
linear relationship between I sp and equivalence ratio was observed between the uniformly 
fueled cases (1, 4 and 8). Very little variation (<1%) in I sp was observed between cases 
of the same overall stoichiometry. These results are summarized in Table 2. 
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Table 2 - Summary of Test Case Performance Results 


Case 

Number 

Overall 

Stoichiometry 

Type of 
Distribution 

Isp 

(s) 

1 

Stoichiometric 

Uniform 

4450 

2 

Stoichiometric 

Transverse 

4429 

3 

Stoichiometric 

Axial 

4459 

4 

Fuel Lean 

Uniform 

4706 

5 

Fuel Lean 

Transverse 

4699 

6 

Fuel Lean 

Axial 

4722 

7 

Fuel Lean 

Partial Fill (Air) 

4676 

8 

Fuel Rich 

Uniform 

4177 

9 

Fuel Rich 

Transverse 

4168 

10 

Fuel Rich 

Axial 

4175 

11 

Fuel Rich 

Partial Fill (H 2 ) 

4204 


Transversely and axially varying fuel/air mixtures remained basically stratified 
throughout the detonation propagation and blowdown process. 


6.2 Conclusions 

The primary conclusion to be drawn from this study is that a lack of fuel-air 
mixing up to a moderate level in a hydrogen fueled airbreathing detonation tube has 
almost no effect on the thrust performance of the system. Put another way, thrust 
performance is almost completely independent of combustion efficiency in the vicinity of 
an overall equivalence ratio of one. This is an encouraging result for the design of such 
systems, as it appears to be unnecessary to go to extremes to achieve good mixing in the 
fuel injection process. These results also imply that it is probably unnecessary to match 
air and fuel valve opening and closing profiles to maintain constant equivalence ratio. 
These results are limited to those cases where the entire fuel/air charge, while not fully 
mixed, is nonetheless fully detonable. 
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There are a number of secondary conclusions that may also be drawn for H 2 /air 
systems. 

1 . Decreased equivalence ratio and partial-filling of the detonation tube are 
essentially equivalent. If anything, it appears slightly more advantageous to lower 
the equivalence ratio as long as the mixture remains detonable. Obviously, a 
point is reached where it is impossible to lower the equivalence ratio and maintain 
a detonable mixture. At this point, further throttling can be achieved by partially 
filling the detonation tube with a lean mixture. 

2. Non-transversely varying mixtures may be treated as one-dimensional, ignoring 
wall shear losses. This means that future studies with such mixtures may be 
safely conducted using one-dimensional codes without fear of significant loss of 
accuracy. This will save considerable time in computational analysis. 

3. Transversely and axially varying mixtures remain essentially stratified throughout 
the detonation propagation and blowdown process. One benefit of this 
observation is that it will be possible to qualitatively evaluate mixing within a 
detonation tube by looking at the distribution of products coming out of the tube, 
thus reducing the need for costly optical test sections in detonation tubes. 

4. H 2 /air detonation tube systems can be adequately computed using simplified 
chemistry, thus reducing required computational time. Early test runs using the 
SPARK code with the full chemical mechanism found in Table A1 showed a 
greater than one order of magnitude increase in computational time required for a 
detonation tube calculation. This was due not only to the increased number of 
species, but also the finer axial grid spacing required to capture with precision the 
trace concentrations of H0 2 and H 2 0 2 in the reactions subsequently eliminated. 


6.3 Recommendations for Future Work 

There are a number of areas worthy of further inquiry that remain at the 
conclusion of this study. These include, 

Increasing the range of equivalence ratios studied. It would be of great value to study 
levels of fuel/air mixture variation that included marginally-detonable and non-detonable 
mixtures. Such mixtures would result from detonation tube film cooling schemes or poor 
mixing configurations. 

1. Performing calculations using a representative hydrocarbon fuel. Hydrocarbon 
fuels are likely candidates for most non-access-to-space applications due to their 
high density and ease of use. It is likely that such systems will behave 
significantly differently. Hydrogen is an excellent monopropellant because of its 
low molecular weight. Therefore, unbumed hydrogen raised to a high 
temperature by the leading shock wave of the detonation front still provides 
significant impulse. A hydrocarbon fuel will not behave in the same manner and 
will likely penalize the system to a greater extent. 
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2. Further development of PDEs is likely to bring forward a number of injection 
patterns not addressed by this study. As realistic patterns become available, 
further analysis should be performed. 

3. Simulations including turbulence should be performed to assess the impact of 
turbulence, both before and after the passage of the detonation wave, on the 
overall system performance. 
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Appendix 

Justification of Truncated Kinetic Mechanism 


The standard H 2 -air kinetic mechanism of Jachimowski (1988) is given in 
Table A1 below. For this study, reactions 8 through 18 have been eliminated. These 
reactions are those that contain H0 2 and H 2 0 2 , species that generally govern the low 
temperature ignition process. However, ignition in a detonation wave occurs 
immediately after the moving shock wave, so the combustible mixture at ignition is not at 
a low temperature. For an equilibrium stoichiometric H 2 -air detonation wave with initial 
conditions of 1 atmosphere and 300 K, the Gordon and McBride (1994) CEA code 
predicts a detonation wave Mach number of 4.8. This corresponds to a temperature and 
pressure behind the shock of approximately 1650 K and 27 atmospheres, respectively, as 
calculated using the CEA code shock option. It is therefore in the neighborhood of this 
temperature and pressure that the necessity of incorporating the H0 2 and H 2 0 2 reactions 
should be considered. 


Table A1 - 18- Step H 2 -Air Mechanism 


Reaction 

Ref.# 

Reaction 

A 

n 

Ea 

(cal/g-mole) 

1 

H 2 +0 2 <=> 2 OH 

0.1700el4 

0.00 

48150 

2 

h + o 2 ^oh + o 

0.1420el5 

0.00 

16400 

3 

h 2 +oh ^ h 2 o+h 

0.3160e08 

1.80 

3030 

4 

o+h 2 <±>oh + h 

0.2070el5 

0.00 

13750 

5 

OH + OH d H 2 0 + 0 

0.5500el4 

0.00 

7000 

6 

OH + H + M d H 2 0 + M 

0.2210e23 

-2.00 

0 

7 

H + H + M <=> H 2 +M 

0.6530el8 

-1.00 

0 

8 

h+o 2 + M « ho 2 +M 

0.3200el9 

-1.00 

0 

9 

oh+ho 2 <=> h 2 o+o 2 

0.5000el4 

0.00 

1000 

10 

H + HO 2 H 2 + 

0.2530el4 

0.00 

700 

11 

H + H0 2 <=> OH + OH 

0.1990el3 

0.00 

1800 

12 

O + HO 2 0 2 + OH 

0.5000el2 

0.00 

1000 

13 

HO 2 + H0 2 o 0 2 + H 2 0 2 

0.1990el3 

0.00 

0 

14 

h 2 +ho 2 <=> h + h 2 o 2 

0.3010el2 

0.00 

18700 

15 

OH + H 2 0 2 <=> H0 2 + H 2 0 

0.1020el4 

0.00 

1900 

16 

h + h 2 o 2 <=> oh+h 2 o 

0.5000el5 

0.00 

10000 

17 

o+h 2 o 2 <=> oh+ho 2 

0.1990el4 

0.00 

5900 

18 

H 2 0 2 +M dOH + OH + M 

0.1210el8 

0.00 

45500 
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The chain-branching reactions in the H 2 /air system are reactions 2, 3, and 4 in 
Table A1 . Kuo (1986) indicates that reaction 2 is the rate-controlling step in the chain- 
branching explosion, as it is the most endothermic. Looking at reactions 2-4 as a stand- 
alone system wherein H 2 and 0 2 are consumed and H 2 0 and H radicals are produced, we 
see that for each H radical consumed, 3 H radicals are produced, for a net of 2 H radicals. 
Thus the net production rate of H radicals in the chain-branching system is 2 times the 
forward rate constant of reaction 2, or 2k t2 . multiplied by the concentrations of H and 0 2 . 
The balancing reaction that prevents the mixture of H 2 and 0 2 from running away into an 
explosion is reaction 8, wherein H0 2 is produced. H0 2 is relatively non-reactive, thus 
terminating the branching reaction chain. Since reaction 8 contains a 3 rd body, it is 
dependent on the number density of molecules, and thus on pressure as well as 
temperature. Therefore, for a given temperature, the importance of the H0 2 producing 
chain terminating reaction is dependant upon the pressure. If we eliminate H0 2 from the 
chemical mechanism, we must determine whether or not we are significantly influencing 
the ignition process. In order to determine this, we need to calculate the temperature at 
which the H radical consumption reaction just balances the H radical production system 
of reactions. If this temperature is significantly below the post-shock temperature, then 
we can safely ignore H0 2 , and consequently H 2 0 2 that is produced from the H0 2 . 

The equality that must be solved for the critical temperature is 
2 k f2 [H][02\ = k fi \H][02][M] (43) 

where [m]= — , as defined by Kuo (1986), the concentration of molecules available 

R u T 

for collisions, i.e. molar density. 

As mentioned previously, each reaction rate constant is of the form 


-Ea 

k f =A*T n *e RJ (28) 

So we can substitute the appropriate values from Table A1 into equation (43) and 
simplify. 


-16400 

(0.284xl0 15 )*<? rj 


(0.320xl0 19 )^ 



(44) 


Solving for temperature at a pressure of 27 atmospheres, the pressure immediately behind 
the shock wave, we find that the creation and consumption of the H radical equilibrate at 
about 1300 K, 350 degrees K below the expected post-shock temperature. Thus we can 
conclude that the H0 2 reactions should not be a significant factor for this study. 
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